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Abstract 

A new method for the reconstruction of the projected mass distribution of clusters of galaxies 
from the image distortion of background galaxies is discussed. This method is essentially 
equivalent to the one we developed previously, i.e., the noise-filtering method, but has several 
practical advantages: (1) It is much easier to implement; (2) it can be easily applied to wide- 
field images, since the constraints on the number of gridpoints are much weaker than for the 
previous method, and (3) it can be easily generalized to more complicated field geometries, 
such as that of the Wide Field Planetary Camera 2 (WFPC2) onboard HST. We have tested 
the performance of our new inversion method (for which a FORTRAN-77 implementation is 
>• | available from the authors) using simulated data, demonstrating that it fares very favourably. 

in : 
o ■ 

1 Introduction 

oo ! 

Q\ ■ The distortion of the images of background galaxies (Tyson, Valdes & Wenk 1990) by the 
tidal gravitational field of clusters of galaxies can be used to obtain a parameter-free re- 
construction of the surface mass density of the cluster (Kaiser & Squires 1993). Several 
modifications of the original reconstruction method were proposed, e.g., to account for dis- 
tortions which are not weak (Seitz & Schneider 1995; Kaiser 1995), to allow an unbiased 
mass reconstruction on a finite field (Schneider 1995; Kaiser et al. 1995; Bartelmann 1995; 
> ! Bartelmann et al. 1996; Seitz & Schneider 1996, hereafter Paper I; Squires & Kaiser 1996), 
and to account for a broad redshift distribution of the background galaxies (Seitz & Schnei- 
der 1997). In this paper, we shall reconsider the second of the above mentioned effects, 
namely mass reconstructions from data on a finite field. In Paper I we have derived a direct 
mass inversion method which is singled out of the infinitely- many unbiased reconstructions 
by identifying a component of the noise (which is due to the intrinsic ellipticity distribution 
of the sources, the discreteness of galaxy images, and observational effects) as such and 
filtering it out. This noise- filter reconstruction has fared very well in numerical simulations 
carried out to compare various finite-field inversions (Paper I; Squires & Kaiser 1996). 

Here, we shall present a slightly revised version of the noise-filter inversion method, 
which removes some of the technical drawbacks of the original formulation. In particular, 
our new method can be applied to arbitrarily-shaped data fields (which is of great interest 
given the geometry of the WF chips of the WFPC2 on-board HST) and can be used with 
better resolution than the previous formulation. In addition, the numerical encoding of the 
new version is substantially easier and requires much less memory. We shall formulate the 
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inversion problem and its solution in Sect. 2, and present some practical issues in Sect. 3. 
Numerical tests of this method in comparison to other reconstruction methods are presented 
in Sect. 4, and we summarize in Sect. 5 our findings. One application of our new method is 
presented in the mass reconstruction of the cluster MS1358+62 by Hoekstra et al. (1998). 

Just before finalizing this manuscript, Lombardi & Bertin (1998) submitted a paper 
to the astro-ph preprint server. Two results of that paper are particularly relevant for the 
present discussion: They have shown that of all (direct) finite-field mass reconstructions, 
those with vanishing curl in the kernel H - see eq. (6) below - have the smallest rms error; 
requiring that noise-free data yield an exact mass reconstruction, they rederived the inversion 
method of Paper I. Second, they have independently derived our new inversion method, 
eq. (7) below, from a variational principle. 



2 Noise-filtered finite-field mass inversion 

We shall assume for simplicity that all source galaxies can be described as being at the same 
redshift; this is not a necessary assumption (see Seitz & Schneider 1997), but simplifies the 
following treatment considerably. Then, let the mass distribution of the cluster be described 
by the dimensionless surface mass density k(0), and the corresponding deflection potential 
be denoted by ip(6), such that the two-dimensional Poisson equation V 2 ip = 2k is satisfied. 
The two components of the complex shear 7 = 71 + i72 are given in terms of the deflection 
potential by 

71 = ^(^,11-^,22), 72 = ^,12, (1) 
where indices separated by a comma denote partial derivatives. 1 The complex reduced shear 

is the expectation value of the observed image ellipticities e, so that the observed image 
ellipticities provide an unbiased estimate of the local value of g, as long as the cluster is 
non-critical. We shall make this assumption here, although it also can be dropped (see 
Seitz & Schneider 1997). As pointed out by Schneider & Seitz (1995), the mass-sheet 
degeneracy (Gorenstein, Falco & Shapiro 1988) allows one to determine (1 — k) only up 
to a multiplicative constant, if no magnification information is used (Broadhurst, Taylor & 
Peacock 1995; Bartelmann & Narayan 1995). Defining 

K(0):=\n[l-n(0)} , (3) 

then K can only be determined up to an additive constant. Kaiser (1995) derived a relation 
between the gradient of K and combinations of first derivatives of g, 



Note that we have changed the sign convention compared to Paper I. 



2 



WK = 



1 - 



-1 



9 



2 



-92 

-92 I + 9i 



( 



91,1+92,2 
92,1 - 91,2 



11(0). 



(4) 



The right-hand-side of this equation can be considered as an observable, obtained from local 
averages of image ellipticities and by finite differencing the resulting field g. 

Equation (4) can be solved (up to an additive constant) by line integration, and several 
schemes for this have been proposed (Schneider 1995; Kaiser et al. 1995; Bartelmann 1995; 
Squires & Kaiser 1996). The reason why different schemes yield different results can be 
seen by noting that the vector field u comes from (noisy) observational estimates, and thus 
will in general not be a gradient field. Therefore, the equation VK = u has no solution in 
general, since u has a rotational component due to observational noise. On the other hand, 
if u is a gradient field, then all line integration schemes are equivalent. 

In Paper I, we split the vector field into a gradient part and a rotational part, 



where s(0) is a scalar field. This decomposition is not unique. However, since the rotational 
component is due solely to noise, we can specify the decomposition uniquely by requiring 
that the mean of rot s over the finite data field U vanishes, and that rot s vanishes if u is 
a gradient field. These two conditions are satisfied if we set s = const on the boundary dll 
of the data field U. Then, identifying VK with VK, the solution of (4) with the rotational 
component removed from u becomes 



where H(0', 0) is a vector field which can be obtained from the Greens function of a Laplace 
equation with Neumann boundary conditions. In Paper I we have derived this equation and 
presented explicit solutions for the cases that U is a circle or a rectangle; in these cases, the 
Greens function can be obtained analytically using geometrical methods. 

Whereas the method of Paper I passed all numerical tests, it has a few features which 
are unwanted: (1) If the geometry deviates from that of a circle or a rectangle, the Greens 
function can no longer be obtained analytically. However, a numerical determination of the 
Greens function is impractical owing to its singularity. For this reason, the mass reconstruc- 
tion of the cluster CI 0939+4713 from WFPC2 data (Seitz et al. 1996) was carried out by 
splitting the field into two rectangles and combine them appropriately in the overlap region. 
This is certainly not the optimal method, since each of the two individual reconstructions 
made no use of the shear information outside the respective rectangle. (2) If a quadratic 
field is covered by an N x N grid of and 0' values, the necessary memory for storing H 
consists of 2N 4 real numbers. Hence, if one increases N beyond ~ 50, the memory require- 
ment quickly approaches the capacity of commonly used workstations. However, due to the 
singularity of the Greens function, one likes to have small grid spacings to obtain an accurate 
numerical estimate of the integral (6) - see Squires & Kaiser (1996) for comments on this 
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point. (3) Although the solution for H was given explicitly in Paper I, it is complicated 
and not easily encoded (though quickly evaluated). In order for the noise-filtering method 
to become a standard and readily available tool, one would like to have an easier method to 
solve for K. 

These three points can be avoided in the following simple manner: Taking the divergence 
of (5) leads to 

V 2 K = V • u . (7a) 

Since s = const on the boundary, rot s is perpendicular to the normal vector n at the 
boundary of U, so that 

n ■ VK = n • u on dU . (76) 

Hence, K can be obtained as the solution of the Neumann problem given by (7a,b). There 
are efficient and quick methods for a numerical solution of this problem; we have employed 
a relaxation method with successive overrelaxation (see Press et al. 1992, p. 857). Choosing 
the overrelaxation parameter as in equation (19.5.21) of Press et al. (1992), a stable solution 
was found after about 20iV iterations on an N x N grid in 6. 

The previously mentioned drawbacks of the method presented in Paper I are now 
avoided. The Neumann problem (7) can be solved for any geometry; for example, for 
the WF-part of the WFPC2 one merely needs to formulate the boundary condition (7b) at 
6 sides, instead of 4 for a rectangle. The memory requirement is reduced to a few times 
iV 2 real numbers, so that N can easily be of order a few hundred. In fact, for N = 200, 
the solution of (7) takes about 2 minutes on an IBM rise 6000 processor. And finally, the 
numerical code for solving (7) shrinks tremendously compared to that needed to evaluate 
H. 



3 Practical implementation 

In order to obtain a mass reconstruction from galaxy ellipticitities, the following three steps 
are needed: 

(1) The galaxy ellipticities are spatially smoothed to obtain an unbiased estimate of the 
local reduced shear. If is the complex ellipticity of the i-th galaxy at position 0^, and A9 
is a smoothing scale, we calculate g at a position as 

where the weight function w is chosen to be 

; 2 

A6 



w(x) = exp I --^2 ) - exp (-q) (9) 



for x < ^yqAO, and zero otherwise. This choice makes w nearly Gaussian and continuous at 
x = ^fqAd, which can be an essential aspect when the derivatives of the components of g 
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need to be evaluated with finite differencing and q is small. In the following we will choose 
q = 9, where (9) becomes almost equal to the case where the correction term is omitted. 
With (8), the reduced shear g can be calculated on a regular grid in 6. 

(2) The vector field u is obtained from g using (4). Finite differencing is employed, 
with one-sided second-order differentiation rules taken at the boundary dU. A further 
differentiation then yields V • u. 

(3) The Neumann problem (7) is then solved, using the method described above. 2 



4 Tests and simulations 

One might wonder whether the mass reconstruction obtained with the method described 
above yields smooth mass profiles. Our method requires differentiation of (noisy) data, so it 
might be suspected that the resulting mass distribution will be quite noisy compared to the 
results of some of the other finite-field inversion methods. These issues have been discussed 
in some detail by Squires & Kaiser (1996); whereas it is not a priori evident that these 
numerical differentiations are unharmful to the resulting mass reconstruction, the numerical 
simulations these authors have carried have shown, in agreement with Paper I, that the 
noise-filter inversion as presented in Paper I yields the least noisy mass estimates of all 
unbiased direct finite-field mass reconstructions that they have tested. 3 Since the method 
proposed here involves a further differentiation of the data, we have to check whether the 
noise level of the reconstruction is not increased by that. 




Fig. 1. Contour- and surface plot for the function —K(x) = — In [1 — k(x)], where k is the two-dimensional 
surface mass density of the lens. The size of the field is 7. '5 x 7/5, and K is calculated on a 50 x 50 grid; 
the spacing in the contour lines is 0.05 



A Fortran 77 code of these three steps is available from the authors by request, both for a rectangular 
data field, and the WFPC2 geometry. 

Inverse methods, such as the maximum likelihood method (Bartelmann et al. 1996) or the maximum 
probability method (Squires & Kaiser 1996) can yield slightly more accurate mass profiles. 
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Fig. 2. Power spectra of the reconstruction error using the mass distribution of Fig. 1, a galaxy density of 50 
per sq. arcminute and a width in the ellipticty distribution of 0.2. The solid, dotted, dashed-dotted, dashed 
and long-dashed curves denote the power spectra obtained using the methods of KS, Paper I (NF), Schneider 
(1995) and the two versions of the newly developed noise-filtering inversion, NF1 and NF2, respectively. The 
gridsize of the final reconstruction was varied from N = 30 to N = 80. 



For this purpose, we have carried out two sets of simulations. In both cases, galaxies 
were distributed randomly on the data field U, with a density of 50 galaxies per square 
arcminute and an intrinsic ellipticity distribution which is assumed to be a Gaussian of 
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width p = 0.2 (see Paper I). In the first set of simulations, a mass distribution for the 
cluster corresponding to the lens model B in Paper I was assumed (see Fig. 1), and the 
'observed' ellipticities were calculated from the intrinsic ellipticities and the local value of 
the reduced shear caused by the lens model. In the second set of simulations, no lens was 
assumed; owing to the intrinsic ellipticity of the sources, the reduced shear as calculated 
from the 'observed' ellipticity does not vanish identically, and so the reconstructed mass 
profile will be different from zero (this is the kind of simulations carried out in Squires & 
Kaiser 1996). The smoothing length was set to AO = 0.'35. 

Reconstructions for the case with a lens were performed using the following methods: 
the original Kaiser & Squires (1993) reconstruction, generalized to account for non-linear 
effects as described in Paper I; a finite-field reconstruction based on line integration (Schnei- 
der 1995); the noise-filtering method as described in Paper I; and the new noise-filtering 
method as presented here. The reconstructions were analyzed by Fourier-decomposition 
of their difference from the input mass distribution (or, more precisely, the input field of 
K). From 50 different realizations of the galaxy population, the power spectrum of this 
difference was obtained, using the same procedure as in Paper I. In Fig. 2, we have plotted 
these power spectra for reconstructions using 40, 80, 30 and 60 gridpoints. Since the KS- 
reconstruction leads to systematic artefacts for finite fields with non-zero shear field outside 
the data region, its power spectrum (solid line) exceeds that of all exact finite field inver- 
sion techniques (see Paper I for more details). The power spectra for the reconstructions 
using the method of Schneider (1995) and the noise-filtering technique of Paper I (NF) are 
shown as dashed-dotted and dotted curves, respectively. The noise filtering inversion devel- 
oped here was implented in two versions (NF1,NF2); in the first case (short-dashed curves) 
(7a&b) was solved on the same grid as for the other inversion techniques. In the second case 
(long-dashed curves), the solution for K was obtained on a grid two times as dense and K 
was estimated on the sparser grid afterwards. 

It turns out that reconstructions with the new technique developed here are always 
smoother than any of the other methods considered here. This is because V • u is used 
instead of u itself. The operation V • u effectively yields a loss of signal and noise on length 
scales of two grid points. To obtain reconstructions with the same resolution we thus double 
the number of gridpoints in NF2, calculate u and its divergence and K on the dense grid. 
Finally K is calculated on the sparser grid by averaging over 4 gridpoints. The power spectra 
of these NF2 reconstructions (long dashed curves) are always very similar to the original 
NF-reconstruction. Those of the NF1 reconstructions are more similar to reconstructions 
on a sparser grid, where the high frequency power is reduced due to the loss of degrees 
of freedom. Since the dotted and long-dashed curves in Fig. 2 are almost the same, this 
demontrates that the recovery of the signal and the sensitivity to the noise in the NF2 and 
NF-method are identical. 

To compare our results with those of Squires & Kaiser (1996) we also consider the case 
with no shear and surface density in the data field (i.e. the 'noise-only-case'). This approach 
investigates the quality of the reconstruction (the quality of the 'no-mass-detection') for the 
case that there is no mass in the field at all, whereas we have investigated before how good a 
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Fig. 3. The same power spectra as in Fig. 2, but without any lens in the field, i.e. k = = 7 



two-dimensional mass distribution can be recovered. Given that a method which involves a 
lot of smoothing will always fare better in the no-lens case than one which spatially resolves 
noise, it is clear that the no-lens comparison is not the relevant test [the "best" inversion in 
that case is obtained by setting H = in (6)!]. 

The curves in Fig. 3 demonstrate again that the noise properties of NF and NF2 recon- 
structions are almost identical, whereas that of NF1 is different for reasons already discussed. 
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Fig. 4. The power spectra of various mass inversion methods, as explained in the main text 

For a dense grid (N = 80) all noise filtering methods become more and more equal, and the 
short wavelength behavior approaches that of KS (solid line). In any case, the KS method 
is by far the best as long as there is no mass in the field. As we already pointed out in 
Paper I, this is because more (and exact) information is used, namely that the shear is (set 
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equal to) zero outside the data field. The fact that the noise of the KS inversion in Fig. 6 of 
Squires & Kaiser (1996) is slightly larger than that of the NF inversion at small wavelengths 
is due to the fact that in their implementation of the KS algorithm, the shear field was not 
obtained by smoothing the galaxy ellipticities, but the inversion was performed by straight 
summation, which leads to shot noise (Seitz & Schneider 1995). 

Squires & Kaiser (1996) suspected that the increase of noise of the finite field inversion 
comes from the fact that they are more sensitive to noise at the boundary of the data field. 
This point is clarified in Fig. 4. The upper and lower solid curves denote the power spectra 
for the NF and KS method on a 40 x 40 grid. The underlying galaxy distribution and thus 
shear field for each of the individual reconstructions is by construction absolutely the same 
for the NF and KS-case. We then embed the true data field U in a two times as large field 
and distribute additional galaxies with the same density and ellipticity distribution in the 
outer region. The galaxies inside U are unchanged. The shear field is calculated in the large 
field and KS-reconstructions are obtained in the same region. We cut out the surface mass 
density in the original field U and calculate the power spectrum of the reconstruction error 
in the same way as for the other mass reconstructions within U. We point out that in this 
case the shear field within U is not the same as in the above case because now galaxies 
outside the field contribute to the estimate of the shear field within U. This makes the 
shear field statistically smoother inside U. But as can be seen in Fig. 4 (long-dashed-dotted 
curve) the reconstruction error within U is larger than for KS-reconstructions of the small 
field (solid line) - because the data outside U are no longer 'ideal assumptions' (7 = 0) but 
noisy measurements affected by the intrinsic ellipticity distribution of the galaxies. We then 
perform reconstructions on the large field where the shear field is obtained in the same way 
as before, but values on gridpoints within U are substituted by the estimate obtained in the 
small field only. Thus the g field and its noise properties within U are now identical to that 
of the KS and NF reconstructions of the small field. At the same time the transition to the 
shear field outside U becomes less continuous which mimics an artifical increase of noise at 
the boundary of U. The power spectrum obtained from KS-reconstructions of that (7-field 
(dotted curve) is higher than the long-dashed-dotted curve, as expected. 

To obtain a KS-reconstruction where almost no information on data outside U is used, 
we increase the noise outside U by doubling the width of the ellipticity distribution for 
galaxies outside U. The shear field is calculated in the large field and the surface density is 
KS-reconstructed. The power spectrum of the reconstruction error within the small field is 
shown as short-dashed-dotted line - and it is very similar to the power spectrum of the finite 
field NF-reconstruction. One could argue that this large increase is caused mainly by the 
fact that by the averaging procedure (8) the increased noise outside U is partly tranferred 
in W. To show that this is not true we again consider the case where the shear field is 
calculated in the large region as before, but where the values inside U are the same as used 
for the KS- and NF-reconstruction of the small field U. We find that the reconstruction 
error is then only marginally decreased (short-dashed curve). But still one could argue that 
in this case the possibly non-smooth transition from the fir-field inside U to that outside 
could significantly contribute to the noise. Therefore we smoothed that transition on the 
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neighboring gridpoints outside the data field. The corresponding power spectrum (long- 
dashed curve) shows that the smootheness of this transition has only a small effect on the 
noise properties of the reconstruction within U. This comparison demonstrates that the KS- 
reconstruction becomes worse the noisier the data outside U are and that the assumption 
of 7 = outside U is responsible for the high quality of the KS-reconstruction if there is no 
mass in the field. Since this is not the case for most fields currently observed, one is urged 
to use a method which is exact on finite fields (see Squires & Kaiser 1996). 

Finally we apply the new noise filtering to the WFPC-2 geometry. Instead of per- 
forming a power-spectrum analysis, we have calculated the mean-square deviation of the 
reconstructed density field K{9) (shifted such that the mean value of K over the field U 
equals the true one) from the input distribution (see also Fig. 10 in Paper I). We consider 
again two cases, the 'no-lens-case' and that of a mass distribution which was now chosen 
similar to that in the Cluster CI 0939 (see Fig. 5) 




Fig. 5. This mass distribution was used when the rms-error for the NF and NF1 were compared; it was 
chosen to similar to that of the cluster C10939. As in Fig. 1 the contours and surface plot shows —K(x) = 
— In [1 — k(x)]. The grid is 40 x 40 and the field of view is 2.5 arcminutes on a side. 

For both cases the galaxy density (60 per square arcminute) , the width of the elliptic- 
ity distribution (p = 0.2) and the smoothing length (A9 = 0.'3) were chosen equal to the 
values for the weak lensing reconstruction of the cluster C10939 (Seitz et al. 1996). The 
reconstructions were obtained on a 40 x 40 grid. In each of the two cases, two different 
reconstructions are analyzed, one where the reconstruction was performed on a square with 
2.' 5 sidelength, and the other where one quarter of the square was removed. Fig. 6 shows 
the rms deviation for these cases, obtained from 500 realizations for each case. For illus- 
tration, only the WFPC-2 part of the square is shown in the first case. When compared 
to reconstructions on the square, the WFPC-2 reconstructions are just slightly more noisy 
close to the additional boundaries of the field, owing to the smaller number of galaxies from 
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which the shear is obtained there. Note that the increase of noise at the 'inner corner' of 
the WFPC-2 is much smaller than that at the 'outer corners', which is due to the fact that 
at the former, more galaxies fall into the filter scale than in the latter case. The increase of 
the noise at those positions where the mass distribution peaks is due to the lack of spatial 
resolution of the inversion, due to the smoothing applied. In contrast to Paper I we have 
not attempted here to adopt an adaptive smoothing, depending on the lens signal, which 
would yield better resolution near the mass peaks. 
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Fig. 6. ... 
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We have derived a new version of the noise-filtering cluster mass reconstruction algorithm 
originally proposed in Paper I, which is easier to implement, easier to use on large fields 
where the required number of gridpoints can quickly exceed the number possible in using 
the method of Paper I, and which can easily be generalized to more complicated geometries; 
the particularly relevant case of the WFPC-2 geometry was considered explicitly. From 
extensive numerical tests we have shown that the noise properties of this version is basically 
identical to that of the method described in Paper I. In agreement with Fig. 6 of Squires & 
Kaiser (1996), we conclude that the noise-filtering method is the best known direct finite-field 
inversion method. The comparison between the maximum probability method (Squires & 
Kaiser 1996) and the method presented here, carried out on the mosaic of WFPC-2 centered 
on the cluster MS1358-I-62 (Hoekstra et al. 1998), yielded no easily visible difference in 
performance of these two methods. 

We thank Bill Press for a very fruitful discussion which triggered the work pre- 
sented here. This work was supported by the " Sonderforschungsbereich 375-95 fur Astro- 
Teilchenphysik" der Deutschen Forschungsgemeinschaft. 
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